﻿// -------------------------------------------------------------------------------------------------------------------------------
// <copyright file="CrcsSparseSymmetricMatrixExtensions.cs" company="Exul">
//     Copyright (c) Exul. All rights reserved.
// </copyright>
// <author>
//     Exul
// </author>
// -------------------------------------------------------------------------------------------------------------------------------

namespace ExulLibrary.Mathematics.LinearAlgebra.Sequential.DoubleComplexExtensions
{
    #region Usings

    using System;
    using System.Globalization;

    using ExulLibrary.Exceptions;
    using ExulLibrary.Mathematics.LinearAlgebra.Exceptions;

    #endregion

    /// <summary>
    ///     Provides extension methods for <see cref="CrcsSparseSymmetricMatrix{T}"/> with
    ///     <see cref="DoubleComplex"/> type parameter.
    /// </summary>
    public static class CrcsSparseSymmetricMatrixExtensions
    {
        #region Methods

        #region public

        /// <summary>
        ///     Adds the second matrix to the specified elements of the first matrix.
        /// </summary>
        /// <param name="augend">
        ///     The augend matrix.
        /// </param>
        /// <param name="augendRowsIndexes">
        ///     The collection of zero-based indexes of rows of augend elements.
        /// </param>
        /// <param name="augendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of augend elements.
        /// </param>
        /// <param name="addend">
        ///     The addend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="augend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="addend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="augendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="augend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="augendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="augend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="augendRowsIndexes"/> and <paramref name="augendColumnsIndexes"/>
        ///     is not in the portrait of <paramref name="augend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendRowsIndexes"/>.Length, <paramref name="addend"/>.RowsCount).
        ///     The number of columns of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendColumnsIndexes"/>.Length, <paramref name="addend"/>.ColumnsCount).
        /// </remarks>
        public static void Add(
            this CrcsSparseSymmetricMatrix<DoubleComplex> augend,
            int[] augendRowsIndexes,
            int[] augendColumnsIndexes,
            DenseSymmetricMatrix<DoubleComplex> addend)
        {
#if SAFE

            if (augend == null)
            {
                throw new ArgumentNullException("augend");
            }

            if (augendRowsIndexes == null)
            {
                throw new ArgumentNullException("augendRowsIndexes");
            }

            if (augendColumnsIndexes == null)
            {
                throw new ArgumentNullException("augendColumnsIndexes");
            }

            if (addend == null)
            {
                throw new ArgumentNullException("addend");
            }

#endif

            var augendSize = augend.RowsCount;
            var begins = augend.Begins;
            var indexes = augend.Indexes;
            var augendDiagonalValues = augend.DiagonalValues;
            var augendValues = augend.TriangleValues;
            var addendValues = addend.Values;

            var addendElementIndex = 0;

            var sizeToAdd = Math.Min(addend.RowsCount, Math.Min(augendRowsIndexes.Length, augendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToAdd; rowIndex++)
            {
                var rowElementIndex = augendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (augendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "augendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, augendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = augendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (augendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "augendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, augendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        augendDiagonalValues[rowElementIndex] += addendValues[addendElementIndex];
                        addendElementIndex++;
                        continue;
                    }

                    var augendElementIndex = rowBegin;
                    while ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] < columnElementIndex))
                    {
                        augendElementIndex++;
                    }

                    if ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] == columnElementIndex))
                    {
                        augendValues[augendElementIndex] += addendValues[addendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);

#endif
                    }

                    addendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Adds the second matrix multiplied by the factor to the specified elements of the first matrix.
        /// </summary>
        /// <param name="augend">
        ///     The augend matrix.
        /// </param>
        /// <param name="augendRowsIndexes">
        ///     The collection of zero-based indexes of rows of augend elements.
        /// </param>
        /// <param name="augendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of augend elements.
        /// </param>
        /// <param name="addend">
        ///     The addend matrix.
        /// </param>
        /// <param name="addendFactor">
        ///     The factor of the addend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="augend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="addend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="augendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="augend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="augendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="augend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="augendRowsIndexes"/> and <paramref name="augendColumnsIndexes"/>
        ///     is not in the portrait of <paramref name="augend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendRowsIndexes"/>.Length, <paramref name="addend"/>.RowsCount).
        ///     The number of columns of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendColumnsIndexes"/>.Length, <paramref name="addend"/>.ColumnsCount).
        /// </remarks>
        public static void Add(
            this CrcsSparseSymmetricMatrix<DoubleComplex> augend,
            int[] augendRowsIndexes,
            int[] augendColumnsIndexes,
            DenseSymmetricMatrix<DoubleComplex> addend,
            DoubleComplex addendFactor)
        {
#if SAFE

            if (augend == null)
            {
                throw new ArgumentNullException("augend");
            }

            if (augendRowsIndexes == null)
            {
                throw new ArgumentNullException("augendRowsIndexes");
            }

            if (augendColumnsIndexes == null)
            {
                throw new ArgumentNullException("augendColumnsIndexes");
            }

            if (addend == null)
            {
                throw new ArgumentNullException("addend");
            }

#endif

            var augendSize = augend.RowsCount;
            var begins = augend.Begins;
            var indexes = augend.Indexes;
            var augendDiagonalValues = augend.DiagonalValues;
            var augendValues = augend.TriangleValues;
            var addendValues = addend.Values;

            var addendElementIndex = 0;

            var sizeToAdd = Math.Min(addend.RowsCount, Math.Min(augendRowsIndexes.Length, augendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToAdd; rowIndex++)
            {
                var rowElementIndex = augendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (augendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "augendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, augendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = augendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (augendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "augendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, augendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        augendDiagonalValues[rowElementIndex] += addendFactor * addendValues[addendElementIndex];
                        addendElementIndex++;
                        continue;
                    }

                    var augendElementIndex = rowBegin;
                    while ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] < columnElementIndex))
                    {
                        augendElementIndex++;
                    }

                    if ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] == columnElementIndex))
                    {
                        augendValues[augendElementIndex] += addendFactor * addendValues[addendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);

#endif
                    }

                    addendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Multiplies the matrix by the dense vector and stores the result in the another vector.
        /// </summary>
        /// <param name="multiplicand">
        ///     The multiplicand matrix.
        /// </param>
        /// <param name="multiplier">
        ///     The multiplier vector.
        /// </param>
        /// <param name="result">
        ///     The vector to store the result.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="multiplicand"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="multiplier"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="result"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="multiplicand"/> to multiply is determined as
        ///     min(<paramref name="multiplicand"/>.RowsCount, <paramref name="result"/>.Length).
        ///     The number of columns of <paramref name="multiplicand"/> to multiply is determined as
        ///     min(<paramref name="multiplicand"/>.ColumnsCount, <paramref name="multiplier"/>.Length).
        /// </remarks>
        public static void Multiply(
            this CrcsSparseSymmetricMatrix<DoubleComplex> multiplicand,
            DenseVector<DoubleComplex> multiplier,
            DenseVector<DoubleComplex> result)
        {
#if SAFE

            if (multiplicand == null)
            {
                throw new ArgumentNullException("multiplicand");
            }

            if (multiplier == null)
            {
                throw new ArgumentNullException("multiplier");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

#endif

            var begins = multiplicand.Begins;
            var indexes = multiplicand.Indexes;
            var multiplicandValues = multiplicand.TriangleValues;
            var diagonalValues = multiplicand.DiagonalValues;
            var multiplierValues = multiplier.Values;
            var resultValues = result.Values;

            var matrixSizeToMultiply = Math.Min(multiplicand.RowsCount, Math.Min(result.Length, multiplier.Length));
            for (var rowIndex = 0; rowIndex < matrixSizeToMultiply; rowIndex++)
            {
                var rowSum = DoubleComplex.Zero;

                var multiplierElement = multiplierValues[rowIndex];

                var elementIndex = begins[rowIndex];
                var rowEnd = begins[rowIndex + 1];
                while (elementIndex < rowEnd)
                {
                    var columnIndex = indexes[elementIndex];

                    rowSum += multiplicandValues[elementIndex] * multiplierValues[columnIndex];
                    resultValues[columnIndex] += multiplierElement * multiplicandValues[elementIndex];

                    elementIndex++;
                }

                rowSum += multiplierElement * diagonalValues[rowIndex];

                resultValues[rowIndex] = rowSum;
            }
        }

        /// <summary>
        ///     Subtracts the second matrix from the specified elements of the first matrix.
        /// </summary>
        /// <param name="minuend">
        ///     The minuend matrix.
        /// </param>
        /// <param name="minuendRowsIndexes">
        ///     The collection of zero-based indexes of rows of minuend elements.
        /// </param>
        /// <param name="minuendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of minuend elements.
        /// </param>
        /// <param name="subtrahend">
        ///     The subtrahend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="minuend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="subtrahend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="minuendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="minuend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="minuendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="minuend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="minuendRowsIndexes"/> and
        ///     <paramref name="minuendColumnsIndexes"/> is not in the portrait of <paramref name="minuend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendRowsIndexes"/>.Length, <paramref name="subtrahend"/>.RowsCount).
        ///     The number of columns of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendColumnsIndexes"/>.Length, <paramref name="subtrahend"/>.ColumnsCount).
        /// </remarks>
        public static void Subtract(
            this CrcsSparseSymmetricMatrix<DoubleComplex> minuend,
            int[] minuendRowsIndexes,
            int[] minuendColumnsIndexes,
            DenseSymmetricMatrix<DoubleComplex> subtrahend)
        {
#if SAFE

            if (minuend == null)
            {
                throw new ArgumentNullException("minuend");
            }

            if (minuendRowsIndexes == null)
            {
                throw new ArgumentNullException("minuendRowsIndexes");
            }

            if (minuendColumnsIndexes == null)
            {
                throw new ArgumentNullException("minuendColumnsIndexes");
            }

            if (subtrahend == null)
            {
                throw new ArgumentNullException("subtrahend");
            }

#endif

            var minuendSize = minuend.RowsCount;
            var begins = minuend.Begins;
            var indexes = minuend.Indexes;
            var minuendDiagonalValues = minuend.DiagonalValues;
            var minuendValues = minuend.TriangleValues;
            var subtrahendValues = subtrahend.Values;

            var subtrahendElementIndex = 0;

            var sizeToSubtract =
                Math.Min(subtrahend.RowsCount, Math.Min(minuendRowsIndexes.Length, minuendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToSubtract; rowIndex++)
            {
                var rowElementIndex = minuendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (minuendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, minuendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = minuendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (minuendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, minuendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        minuendDiagonalValues[rowElementIndex] -= subtrahendValues[subtrahendElementIndex];
                        subtrahendElementIndex++;
                        continue;
                    }

                    var minuendElementIndex = rowBegin;
                    while ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] < columnElementIndex))
                    {
                        minuendElementIndex++;
                    }

                    if ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] == columnElementIndex))
                    {
                        minuendValues[minuendElementIndex] -= subtrahendValues[subtrahendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);

#endif
                    }

                    subtrahendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Subtracts the second matrix multiplied by the factor from the specified elements of the first matrix.
        /// </summary>
        /// <param name="minuend">
        ///     The  minuend matrix.
        /// </param>
        /// <param name="minuendRowsIndexes">
        ///     The collection of zero-based indexes of rows of minuend elements.
        /// </param>
        /// <param name="minuendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of minuend elements.
        /// </param>
        /// <param name="subtrahend">
        ///     The subtrahend matrix.
        /// </param>
        /// <param name="subtrahendFactor">
        ///     The factor of the subtrahend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="minuend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="subtrahend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="minuendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="minuend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="minuendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="minuend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="minuendRowsIndexes"/> and
        ///     <paramref name="minuendColumnsIndexes"/> is not in the portrait of <paramref name="minuend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendRowsIndexes"/>.Length, <paramref name="subtrahend"/>.RowsCount).
        ///     The number of columns of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendColumnsIndexes"/>.Length, <paramref name="subtrahend"/>.ColumnsCount).
        /// </remarks>
        public static void Subtract(
            this CrcsSparseSymmetricMatrix<DoubleComplex> minuend,
            int[] minuendRowsIndexes,
            int[] minuendColumnsIndexes,
            DenseSymmetricMatrix<DoubleComplex> subtrahend,
            DoubleComplex subtrahendFactor)
        {
#if SAFE

            if (minuend == null)
            {
                throw new ArgumentNullException("minuend");
            }

            if (minuendRowsIndexes == null)
            {
                throw new ArgumentNullException("minuendRowsIndexes");
            }

            if (minuendColumnsIndexes == null)
            {
                throw new ArgumentNullException("minuendColumnsIndexes");
            }

            if (subtrahend == null)
            {
                throw new ArgumentNullException("subtrahend");
            }

#endif

            var minuendSize = minuend.RowsCount;
            var begins = minuend.Begins;
            var indexes = minuend.Indexes;
            var minuendDiagonalValues = minuend.DiagonalValues;
            var minuendValues = minuend.TriangleValues;
            var subtrahendValues = subtrahend.Values;

            var subtrahendElementIndex = 0;

            var sizeToSubtract =
                Math.Min(subtrahend.RowsCount, Math.Min(minuendRowsIndexes.Length, minuendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToSubtract; rowIndex++)
            {
                var rowElementIndex = minuendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (minuendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, minuendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = minuendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (minuendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, minuendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        minuendDiagonalValues[rowElementIndex] -= subtrahendFactor * subtrahendValues[subtrahendElementIndex];
                        subtrahendElementIndex++;
                        continue;
                    }

                    var minuendElementIndex = rowBegin;
                    while ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] < columnElementIndex))
                    {
                        minuendElementIndex++;
                    }

                    if ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] == columnElementIndex))
                    {
                        minuendValues[minuendElementIndex] -= subtrahendFactor * subtrahendValues[subtrahendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);

#endif
                    }

                    subtrahendElementIndex++;
                }
            }
        }

        #endregion

        #endregion
    }
}